home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 20 / Cream of the Crop 20 (Terry Blount) (1996).iso / math / alged34.zip / ALGEDSRC.ZIP / ALGEXPR.C < prev    next >
C/C++ Source or Header  |  1996-06-06  |  20KB  |  787 lines
  1. /*--------------------------------------------------------------------
  2.    Alged:  Algebra Editor    henckel@vnet.ibm.com
  3.  
  4.    Copyright (c) 1994 John Henckel
  5.    Permission to use, copy, modify, distribute and sell this software
  6.    and its documentation for any purpose is hereby granted without fee,
  7.    provided that the above copyright notice appear in all copies.
  8. */
  9. #include "alged.h"
  10. /*--------------------------------------------------------------------
  11.    associate
  12.    rotate the association on add,sub,mul,div
  13.    a+b-c => a-c+b => b+a-c => b-c+a
  14.    note, we make use of the fact that ADD+1 = SUB and ADD is even.
  15. */
  16. void associate(node *p) {
  17.   node *a,*b;
  18.   int opr;
  19.  
  20.   opr = p->kind;
  21.   b = p->rt;
  22.  
  23.   if (opr==DIV && b->kind==MUL) {     /* special handle for bisected */
  24.     p->rt = b->rt;
  25.     b->rt = b->lf;
  26.     b->lf = p->lf;
  27.     p->lf = b;
  28.     b->kind = DIV;
  29.     b = p->rt;
  30.   }
  31.   if (opr==ADD || opr==MUL || opr==SUB || opr==DIV) {
  32.     a = p;
  33.     while ((a->lf->kind|1) == (opr|1)) {
  34.       a->kind = a->lf->kind;
  35.       a->rt = a->lf->rt;
  36.       a = a->lf;
  37.     }
  38.     a->kind = opr;
  39.     if (opr&1) {           /* subtr or divide */
  40.       a->rt = b;
  41.     }
  42.     else {
  43.       a->rt = a->lf;
  44.       a->lf = b;
  45.     }
  46.   }
  47.   else if (opr==EQU) {           /* commute equality */
  48.     p->rt = p->lf;
  49.     p->lf = b;
  50.   }
  51. }
  52.  
  53. /*--------------------------------------------------------------------
  54.    commute
  55.    commute on add,sub,mul,div
  56.    e.g.   a/b ==>  b^(-1)/a^(-1)
  57. */
  58. void commuteNOTUSED(node *p) {
  59.   node *a,*b;
  60.  
  61.   a = p->lf;
  62.   b = p->rt;
  63.   if (p->kind==MUL || p->kind==ADD || p->kind==EQU) {
  64.     p->lf=b; p->rt=a;
  65.   }
  66.   else if (p->kind==DIV || p->kind==SUB) {     /* non-commutative */
  67.     if (a->kind==p->kind+1 &&
  68.         b->kind==p->kind+1 &&
  69.         a->rt->kind==NUM &&
  70.         b->rt->kind==NUM) {              /* has coefficients */
  71.       p->lf=b; p->rt=a;
  72.       a->rt->value = -a->rt->value;
  73.       b->rt->value = -b->rt->value;
  74.     }
  75.     else {                          /* need to make coefficients */
  76.       p->lf = newoper(p->kind+1);
  77.       p->rt = newoper(p->kind+1);
  78.       p->lf->lf = b;
  79.       p->lf->rt = newnum(-1);
  80.       p->rt->lf = a;
  81.       p->rt->rt = newnum(-1);
  82.     }
  83.   }
  84. }
  85.  
  86. /*--------------------------------------------------------------------
  87.       a^x * a^y ==> a^(x+y)
  88.       x^a * y^a ==> (xy)^a
  89. */
  90. int expjoin(node *p) {
  91.   int i,r=0;
  92.   node *a,*b;
  93.  
  94.   for (i=0; i<p->nump; ++i)
  95.     r+=expjoin(p->parm[i]);
  96.  
  97.   a = p->lf;
  98.   b = p->rt;
  99.   if (p->kind==MUL || p->kind==DIV) {
  100.     if (a->kind==EXP && b->kind==EXP &&
  101.         equal(a->lf,b->lf)) {      /*  a^x * a^y = a^(x+y) */
  102.       freetree(a->lf);
  103.       a = a->rt;
  104.       freenode(p->lf);
  105.       b->kind = p->kind-2;
  106.       p->kind = EXP;
  107.       p->lf = b->lf;
  108.       b->lf = a;
  109.     }
  110.     else if (b->kind==EXP &&
  111.         equal(a,b->lf)) {               /*  a * a^y = a^(1+y) */
  112.       freetree(a);
  113.       a = newnum(1);
  114.       b->kind = p->kind-2;
  115.       p->kind = EXP;
  116.       p->lf = b->lf;
  117.       b->lf = a;
  118.     }
  119.     else if (a->kind==EXP &&
  120.              equal(a->lf,b)) {          /*  a^x * a = a^(x+1) */
  121.       freetree(b);
  122.       b = newnum(1);
  123.       a->kind = p->kind-2;
  124.       p->kind = EXP;
  125.       p->lf = a->lf;
  126.       a->lf = a->rt;
  127.       a->rt = b;
  128.       p->rt = a;
  129.     }
  130.     else if (a->kind==EXP && b->kind==EXP &&
  131.         equal(a->rt,b->rt)) {      /*  x^a * y^a = (xy)^a */
  132.       freetree(b->rt);
  133.       p->rt = a->rt;
  134.       a->rt = b->lf;
  135.       freenode(b);
  136.       a->kind = p->kind;
  137.       p->kind = EXP;
  138.     }
  139.     else return r;
  140.     ++r;
  141.   }
  142.   return r;
  143. }
  144.  
  145. /*--------------------------------------------------------------------
  146.    remove sub and div
  147.  
  148.    x/y  =  x*y^-1
  149.    x-y  =  x+y*-1
  150. */
  151. int nosubdiv(node *p) {
  152.   int i,r=0;
  153.   node *a,*b;
  154.  
  155.   for (i=0; i<p->nump; ++i)
  156.     r+=nosubdiv(p->parm[i]);
  157.  
  158.   a = p->lf;
  159.   b = p->rt;
  160.   if (p->kind==SUB || p->kind==DIV) {
  161.     ++r;
  162.     --p->kind;
  163.     if (b->kind==NUM && p->kind==ADD)           /* a-2 = a+(-2) */
  164.        b->value = -b->value;
  165.     else if (b->kind==NUM && !whole(b->value))
  166.        b->value = 1.0/b->value;
  167.     else {
  168.       a = newoper(p->kind+2);    /* MUL or EXP */
  169.       a->lf = b;
  170.       a->rt = newnum(-1);
  171.       p->rt = a;
  172.     }
  173.   }
  174.   return r;
  175. }
  176.  
  177. /*--------------------------------------------------------------------
  178.    bisect node
  179.  
  180.    This converts expressions to canonical form in which
  181.    1. + * are left assoc, ^ is right assoc.
  182.    2. scope of / is increased, e.g. (a/b)*c ==> ac/b,
  183.    3. x+y*-2 is changed to x-y*2
  184.    3' x+(-2) is changed to x-2
  185.    4. x*y^-2 is changed to x/y^2
  186.    4' x*(0.5) is changed to x/2
  187.    5. scope of ^ is reduced
  188.    The name "bisect" means that each MUL clag is broken into two pieces the
  189.    numerator elements and the denominator ones.
  190. */
  191. int bisect(node *p) {
  192.   int i,r=0;
  193.   node *a,*b;
  194.  
  195.   for (i=0; i<p->nump; ++i)
  196.     r+=bisect(p->parm[i]);
  197.  
  198.   a = p->lf;
  199.   b = p->rt;
  200.   switch (p->kind) {
  201.   case ADD:
  202.   case MUL:
  203.     /*--------------------------------------------------------------------
  204.        Add or Multiply
  205.     */
  206.     if (b->kind==p->kind ||              /* a+(b+c) = a+b+c */
  207.         b->kind==p->kind+1) {            /* a+(b-c) = a+b-c */
  208.       swingb;
  209.       i = b->kind;
  210.       b->kind = p->kind;
  211.       p->kind = i;
  212.     }
  213.     else if (p->kind==ADD &&             /* a+(-2) = a-2 */
  214.              b->kind==NUM &&
  215.              b->value < 0) {
  216.       p->kind = SUB;
  217.       b->value = -b->value;
  218.     }
  219.     else if (b->kind==p->kind+2 &&       /* a+b*-2 = a-b*2 */
  220.         b->rt->kind==NUM &&
  221.         b->rt->value < 0) {
  222.       ++p->kind;
  223.       b->rt->value = -b->rt->value;
  224.     }
  225.     /*--------------------------------------------------------------------
  226.        Add only
  227.     */
  228.     else if (p->kind==ADD &&
  229.         aop(b->kind)) {                   /* a+(b-c) = (a+b)-c */
  230.       swingb;
  231.       p->kind = b->kind;
  232.       b->kind = ADD;
  233.     }
  234.     /*--------------------------------------------------------------------
  235.        Multiply only
  236.     */
  237.     else if (p->kind==MUL &&
  238.         a->kind==p->kind+1) {       /* a-b+c = a+c-b */
  239.       p->rt = a->rt;
  240.       a->rt = b;
  241.       ++p->kind;
  242.       --a->kind;
  243.     }
  244.     else if (p->kind==MUL &&
  245.         a->kind==p->kind+2 &&       /* b*-2+a = a-b*2 */
  246.         a->rt->kind==NUM &&
  247.         a->rt->value < 0) {
  248.       ++p->kind;
  249.       p->lf = b;
  250.       p->rt = a;
  251.       a->rt->value = -a->rt->value;
  252.     }
  253.     else break;
  254.     ++r; break;
  255.   case SUB:
  256.   case DIV:
  257.     /*--------------------------------------------------------------------
  258.        Subtract or divide
  259.     */
  260.     if (b->kind==p->kind+1 &&       /* a-b*-2 = a+b*2   NOT NECES. */
  261.         b->rt->kind==NUM &&
  262.         b->rt->value < 0) {
  263.       --p->kind;
  264.       b->rt->value = -b->rt->value;
  265.     }
  266.     /*--------------------------------------------------------------------
  267.        Subtract only
  268.     */
  269.     else if (p->kind==SUB &&             /* a-(-2) = a+2 */
  270.         b->kind==NUM &&
  271.         b->value < 0) {
  272.       p->kind = ADD;
  273.       b->value = -b->value;
  274.     }
  275.     else if (p->kind==SUB &&
  276.          aop(b->kind)) {                   /* a-(b+c) = (a-b)-c */
  277.       swingb;
  278.       p->kind = (ADD+SUB) - b->kind;
  279.       b->kind = SUB;
  280.     }
  281.     /*--------------------------------------------------------------------
  282.        Divide only
  283.     */
  284.     else if (p->kind==DIV &&
  285.         b->kind==p->kind) {             /* a-(b-c) = a+c-b */
  286.       p->lf = b;
  287.       p->rt = b->lf;
  288.       b->lf = a;
  289.       --b->kind;
  290.     }
  291.     else if (p->kind==DIV &&
  292.         a->kind==p->kind) {        /* a-b-c = a-(b+c) */
  293.       swinga;
  294.       --a->kind;
  295.     }
  296.     else break;
  297.     ++r; break;
  298.   case EXP:
  299.     /*--------------------------------------------------------------------
  300.        exponent
  301.     */
  302.     if (a->kind==EXP) {                  /* (x^y)^z = x^(y*z) */
  303.       swinga;
  304.       a->kind = MUL;
  305.     }
  306.     else if (b->kind==NUM &&             /* a^(-2) = 1/a^2 */
  307.              b->value < 0) {
  308.       p->kind = DIV;
  309.       p->lf = newnum(1);
  310.       p->r